function [model,errorMes] = solveQTSM_constRsmooth(alpha,rhor,beta,psi,mu,phi,sigma,setup)

%% Allocating memory 
maxMat     = max(setup.matSelect);            % max maturity
ny         = size(setup.matSelect,2);         % number of bond yields
nxtil      = 2*setup.nm + setup.nn;           % number of pricing factors

% Law of motion for factors under Q 
h0Q = phi*mu;
hxQ = eye(nxtil)-phi;

% Check if co-variance matrix for innovations is positive semi-definite
sigma2 = sigma*sigma';
if any(eig(sigma2)<0)
    errorMes = 1;
    model = [];
    return
end
%% Compute bond pricing coefficients
% Initialise memory for bond price coefficients in state 1
A1 = zeros(1,maxMat);
B1 = zeros(nxtil,maxMat);
C1 = zeros(nxtil,nxtil,maxMat);
D1 = zeros(1,maxMat);

% Compute recursive coefficients for bond prices
for n = 1:maxMat   
    % One-period bond prices use the initial conditions
    if n == 1
        % Compute the recursive bond pricing coefficients.  
        A1(1,n)   = - alpha;
        B1(:,n)   = - beta;
        C1(:,:,n) = - psi;
        D1(1,n)   = - rhor;

    % Longer bond prices depend on the previous maturity
    else
        gamma    = eye(nxtil)-2*sigma'*C1(:,:,n-1)*sigma;
        %invgamma = inv(gamma);
        invgamma = gamma\eye(nxtil);
        detGamma = det(gamma);
        if detGamma < 0
            errorMes = 1;
            model = struct();
            return
        end
        A1(1,n) = A1(1,n-1)+ (D1(1,n-1)-1)*alpha + B1(:,n-1)'*h0Q ...
            + h0Q'*C1(:,:,n-1)*h0Q ...
            + 0.5*B1(:,n-1)'*sigma*invgamma*sigma'*B1(:,n-1) ...
            + 2*h0Q'*C1(:,:,n-1)*sigma*invgamma*sigma'*B1(:,n-1)...
            + 2*h0Q'*C1(:,:,n-1)*sigma*invgamma*sigma'*C1(:,:,n-1)'*h0Q ...
            - 0.5*log(detGamma);       
        
        B1(:,n) = (B1(:,n-1)'*(eye(nxtil)-phi)+(D1(1,n-1)-1)*beta'...
                    +2*h0Q'*C1(:,:,n-1)*hxQ ...
                    +2*B1(:,n-1)'*sigma*invgamma*sigma'*C1(:,:,n-1)'*hxQ...
                    +4*h0Q'*C1(:,:,n-1)*sigma*invgamma*sigma'*C1(:,:,n-1)'*hxQ)';     

        C1(:,:,n) = hxQ'*C1(:,:,n-1)*hxQ + (D1(1,n-1)-1)*psi ...
                    + 2*hxQ*C1(:,:,n-1)*sigma*invgamma*sigma'*C1(:,:,n-1)'*hxQ;
                
        D1(1,n)   = (D1(1,n-1)-1)*rhor;

    end
end


%% Reporting output in the desired format: annualized yields in pct.
% Spot: y_{n,t} = -(1/n)*p_{n,t}
nx         = nxtil;

% The short rate
r0   = 12*(-A1(1)/1);
rx   = 12*(-B1(:,1)'/1);
rxx  = 12*(-reshape(C1(:,:,1),1,nx^2)/1);
g0   = 12*(-A1(1,setup.matSelect)'./setup.matSelect');
gx   = 12*(-B1(:,setup.matSelect)'./repmat(setup.matSelect',1,nx));
gxx  = zeros(ny,nx^2);
gr   = 12*(-D1(1,setup.matSelect)'./setup.matSelect');
for i=1:ny
    gxx(i,:) = 12*(-reshape(C1(:,:,setup.matSelect(1,i)),nx^2,1)./repmat(setup.matSelect(1,i),1,nx^2)');
end

% Storing output
model = struct('ny',length(g0),'nx',nx,'g0',g0,'gx',gx,'gxx',gxx,'h0Q',h0Q,'hxQ',hxQ,'matSelect',setup.matSelect,...
    'r0',r0,'rx',rx,'rxx',rxx,...
    'alpha',alpha,'beta',beta,'psi',psi,'phi',phi,'mu',mu,'sigmaxtil',sigma,...
    'A',A1,'B',B1,'C',C1);

% Report no error if applicable
errorMes = 0;


end